####This script will calculate average reptile abundance and standard deviation on a site by site basis for 4 species of skink

#Clear workspace

rm(list=ls())

in.dir = "C:/R/In/"
out.dir = "C:/R/Out/"

#Set working directory
setwd(in.dir)

#load necessary libraries

#list the libraries needed
necessary=c("adehabitat","SDMTools","raster")                                
#check if library is installed
installed = necessary %in% installed.packages()
#if library is not installed, install it
if (length(necessary[!installed]) >=1) install.packages(necessary[!installed], dep = T)

#load the libraries
for (lib in necessary) library(lib,character.only=T)

#Import .csv file containing the number of individual skinks counted per sample and save as object
from.mvd = read.csv("Master Reptile Abundance Dataset.csv", header = TRUE)

# First adjust raw data so that number of animals per survey is changed to reflect abundance in individuals per hectare (NB: One survey is presumed to cover .785 ha)

reptile = data.frame(georef_ID=from.mvd$georef_ID, east=from.mvd$east, north=from.mvd$north, LAMCOGG=(from.mvd$LAMCOGG*.785), CARRUBR=(from.mvd$CARRUBR*.785), SAPBASI=(from.mvd$SAPBASI*.785), GNYQUEE=(from.mvd$GNYQUEE*.785))

#Create a function to calculate the Mean Absolute Deviation of a distribution

t.mean = function(x){return(mean(x,na.rm=T))}
t.sd = function(x){return(sd(x,na.rm=T))}

#Aggregate reptile datasets to perform calculation of mean and SD skink abundance for multiple surveys at individual sites

LAMCOGGmean = aggregate(reptile$LAMCOGG, by=list(georef_ID=reptile$georef_ID), FUN = t.mean)
LAMCOGGsd= aggregate(reptile$LAMCOGG, by=list(georef_ID=reptile$georef_ID), FUN = t.sd)
CARRUBRmean = aggregate(reptile$CARRUBR, by=list(georef_ID=reptile$georef_ID), FUN = t.mean)
CARRUBRsd= aggregate(reptile$CARRUBR, by=list(georef_ID=reptile$georef_ID), FUN = t.sd)
SAPBASImean = aggregate(reptile$SAPBASI, by=list(georef_ID=reptile$georef_ID), FUN = t.mean)
SAPBASIsd= aggregate(reptile$SAPBASI, by=list(georef_ID=reptile$georef_ID), FUN = t.sd)
GNYQUEEmean = aggregate(reptile$GNYQUEE, by=list(georef_ID=reptile$georef_ID), FUN = t.mean)
GNYQUEEsd= aggregate(reptile$GNYQUEE, by=list(georef_ID=reptile$georef_ID), FUN = t.sd)

#Create a dataframe that is georef_ID,east,north,and each of the 4 measures of abundance & SD

master = data.frame(georef_ID=LAMCOGGmean$georef_ID, LAMCOGGmean$x, LAMCOGGsd$x, CARRUBRmean$x, CARRUBRsd$x, SAPBASImean$x, SAPBASIsd$x, GNYQUEEmean$x, GNYQUEEsd$x)

#Export dataframe to import back into MVD and line up georef_ID with east & north.  Then re-import file and intersect with Environmental Suitability Values at each site for each species

write.csv(x=master,file=paste(out.dir, "abundancebysite.csv", sep=""), row.names=F)

#Reimport 'master' with east/north from MVD

abund = read.csv("Abundance By Site.csv", header=TRUE)

#Intersect 'abund' with environmental suitability values from maxent outputs, first create a paired list of points to intersect.
#In this case, no loop was written, the name of the ASCII file and column to append to object 'pnts' was manually changed for each of 4 species CARRUBR, LAMCOGG, SAPBASI, and GNYQUEE
t.file = paste(in.dir, "ASCII/","LAMCOGG.asc", sep ="")
t.asc = read.asc(t.file)
east = c(abund$east)
north = c(abund$north)
pnts = data.frame(east=east, north=north)
pnts$LAMCOGGmean = join.asc(pnts, t.asc)
rm(t.asc)
final = merge(pnts,abund,by=c("east","north"))

#Write out a new data frame with new column names

final2 = data.frame(georef_ID=final$georef_ID, east=final$east, north=final$north, carrubr_abund=final$CARRUBRmeanx, carrubr_sd=final$CARRUBRsdx, carrubr_es=final$CARRUBRmean, lamcogg_abund=final$LAMCOGGmeanx, lamcogg_sd=final$LAMCOGGsdx, lamcogg_es=final$LAMCOGGmean, gnyquee_abund=final$GNYQUEEmeanx, gnyquee_sd=final$GNYQUEEsdx, gnyquee_es=final$GNYQUEEmean, sapbasi_abund=final$SAPBASImeanx, sapbasi_sd=final$SAPBASIsdx, sapbasi_es=final$SAPBASImean)

#Write out the final data frame as a .csv file

write.csv(x=final2,file=paste(out.dir,"pairedabundanceandenvironmentalsuitabilitybysite.csv",sep=""), row.names=F)
